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The paper describes two general problems encountered in computational assignments at the in- 
troductory level. First, novice students often treat computer code as almost magic incantations, 
and like novices in many fields, have trouble creating new algorithms or procedures to solve novel 
problems. Second, the nature of computational studies often means that the results generated are 
interpreted via theoretically devised quantities, which may not meet a student's internal standards 
for proof when compared to an experimental measurement. 

The paper then offers a lab/programming assignment, used in a calculus-based physics course, 
which was devised to address these problems. In the assignment, students created a computational 
model of the thermal energy transfer involved in heating an iron rod in a blacksmith's forge. After 
creating the simulation, students attended a blacksmithing seminar and had a chance to work with 
iron and take data on its heating in a coke-fueled forge. On their return to campus, students revised 
their computational models in light of their experimental data. 



I. INTRODUCTION 

Computational Physics is quickly becoming a standard 
component of many departments' degree program and 
the utility of the field cannot be denied. It's hard to imag- 
ine a modern physics lab without a computer to perform 
data analysis, and more than just making the task more 
bearable, most present experimental programs would be 
impossible without computational automation. [1, 2] 

In addition, the ability to create a working model 
of a "real" system by distilling out the essential physi- 
cal/mathematical relationships, and then implementing 
these relationships in a way that captures part of na- 
ture is an existential thrill for many students [3, 4, 5]. 
When we allow our students to tackle novel problems 
and discover real truths on their own, via a simulation 
they create, we recruit new physics majors [6, 7]. 

A. Novice and Expert Behavior 

When learning how to program, students often seem to 
grab hold of a certain example program and then treat 
the syntax therein as a magic incantation. So long as they 
don't deviate too far from the example code, their pro- 
gram should compile, run, and produce something close 
to the desired result. In so doing, students miss the ac- 
tual structure of the computer language and the ability 
to have a computer do "whatever you want" rather than 
"only the things I have examples for" [8]. 

This is not a trait specific to computational physics, or 
even computer programming in general. Most novices in 
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any field start the learning process by copying the actions 
of an expert. In physics, this often occurs in homework 
problems where students are sometimes asked to parrot 
the professor's class examples in different contexts. In 
the same vein, an infant of about six month's age, when 
seated at the dinner table, will start opening and closing 
their mouth, just like their parents do when they lift a 
utensil to their mouth. 

A more sophisticated student will make a plan for how 
the program should look and function before actually 
writing the code. Instead of looking through their li- 
brary of example programs for one that does something 
close to the job at hand, an expert will imagine the al- 
gorithmic steps necessary to solve the problem, and then 
figure out what computational machinery (variables, op- 
erations, and syntax) best suit the algorithm. [9] 

This is not a new division of skill. The novice in in- 
troductory physics looks at a test problem about the 
flight of a baseball and tries to figure out what worked 
class example best matches the problem, perhaps pick- 
ing the worked trebuchet example from class. An expert, 
working the same problem, will likely exhibit clustering 
or "chunking" of knowledge [11]. Rather than thinking 
about specific worked solutions, experts tend to think of 
all the pieces of physics related to the question. In this 
case, the expert might call to mind kinematics, energy 
conservation, air resistance, and Newton's second law. 
Although the expert's solution might take longer to pro- 
duce (and contain extraneous details), the fact that the 
problem is novel is not a insurmountable difficulty for 
an expert, [10]. By contrast, a novice can generally only 
solve problems that are in their library of solutions. 

The "best ways" to move students to the expert state 
arc lcgion[12, 13, 14], and there's no utility in repeating 
the list here. In the context of programming though, a 
good start towards creating expert behavior is to require 
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it [15, 16, 17]. If we assume that the instruction and 
materials are clear and sufficient, then the first law of 
teaching, "If you want students to do something, grade 
it!" [18] would seem to apply. 



B. Misplaced Trust in Technology 

Most calculus textbooks contain a section in the first 
chapter that might as well be called "Lies my calculator 
told me" [19]. As most instructors well know, novice stu- 
dents often trust calculators more than their own heads, 
and this trust is dangerous, particularly because the cal- 
culator is only as reliable as the person who designed and 
programmed it. Of course, very few physicists see value 
in requiring students to be able to compute something 
like log(1.9723) by hand, rather, the idea is that stu- 
dents should have the habit of thinking critically about 
the answer a calculator produces. 

As a more concrete example, consider the student who 
is using a calculator to evaluate the angle 9 produced 
by the opposite and hypotenuse sides of a triangle. If 
the hypotenuse is h = 1.0, and the opposite side is b = 
0.2, 0.4, 0.6, 0.8, ... the student will have little trouble with 
the arcsin button on their calculator (using the relation 
9 = arcsin £). If however, the student notices a trend in 
the value of 9, and out of curiosity plugs in the value b = 
1.2, their calculator will fail them. Of course, a very good 
calculator might map arcsin to the complex plane and 
produce the answer 9 = 1.57— 0.62z, but most calculators 
would just produce "ERROR". A novice student will 
generally try the computation again, get the same result, 
and then proclaim that their calculator is "broken" . A 
more sophisticated student would hopefully think about 
the geometry implied by b = 1.2, and perhaps generalize 
that for values of b > 1.0 the triangle is no longer right 
and will return meaningless results for 9. 

Speaking generally, expert students think critically 
about the answers they generate, either by hand, by cal- 
culator, or by computer, and use unexpected answers ei- 
ther as an indication of error, or as an indicator of some 
deeper sophistication to the problem [16, 17]. In compu- 
tational physics this trait is a complicated one to rein- 
force, because the nature of the systems studied generally 
means that the systems are too small, large, hot, cold, 
or expensive to watch experimentally [20]. The common 
approach is to monitor some other quantity, like the sys- 
tem's net momentum or total energy, and check to see 
if the simulated behavior matches with the theoretically 
postulated result. 

For example, in an idealized simulation of the moon's 
orbit around the earth, a student can monitor the grav- 
itational potential energy of the moon and the kinetic 
energy of the moon, and compare this observation with a 
theoretical conservation law. A similar measurement of 
the net momentum of the moon would reveal that this 
quantity is not conserved, which an advanced student 
might link to the external gravitational force exerted on 



the earth by the moon. 

The problem with such measurements is that they're 
the result of theoretical proofs or conjectures, rather 
than experimental observations. In the case of most 
introductory-level students I have worked with, a plot 
of kinetic energy is not as persuasive as a computer gen- 
erated temperature that corresponds to an actual ther- 
mometer measurement, taken in person. Conservation 
laws are beautiful, but our beginning students don't be- 
lieve them. Many college freshmen are still in the con- 
crete Piagetian stage [21], and accordingly, "proof" to an 
introductory physics student needs to include something 
a student can touch. 



II. GOAL OF THE PRESENT WORK 

To summarize, in computational physics an ideal prob- 
lem will have solution which can be compared to "real" 
data which ideally, is collected by the students. It seems 
that agreement between student taken data and a stu- 
dent written simulation is more powerful for learning that 
agreement with "facts from a book." An ideal problem 
will also avoid the problem of "solution by finding a suit- 
able example," described above, by requiring a student to 
translate physics and mathematics into computer code, 
which forces a student to think about what the mathe- 
matics and numerical approximation actually mean (and 
thus moves a student along the path toward expert think- 
ing). 

The rest of this paper describes a computational exer- 
cise which is intended to address the two problems de- 
scribed above. The problem involves a student- written 
simulation of heat flow within a bar of iron, which the 
class then compared to an experiment in the field at a 
blacksmithing workshop. After returning from the work- 
shop, students modified their simulations to better ac- 
count for the experimental data they recorded. This 
project was given in the second semester of a Univer- 
sity Physics course at Winona State University. The class 
was populated mainly with engineering students who had 
been exposed to programming via the VPython exercises 
in the "Matter and Interactions" [22, 23] introductory 
physics text. Rather than simply using a computer to 
determine a numerical answer, this class had an empha- 
sis on trying to learn about the world through computer 
simulation. To make this goal more than just a slogan, 
computer programming assignments in this class were 
evaluated with the rubric in appendix A. 



III. THE BLACKSMITHING PROBLEM 

The specific problem used in class is described in the 
following section. The general idea was for students to 
get more out of a field experience by taking the time 
to think about and write a numerical model before go- 
ing on the trip. After returning, the students had a 
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chance to reflect on the validity of their numerical work 
by comparing computational predictions to actual exper- 
imental measurements. This approach of prediction — > 
data collection — > reflection is not novel, [24, 25], but 
real comparisons of simulation and experimental data 
seem somewhat rare in introductory computational sci- 
ence. 



A. Background and Preparation 

With the thought that other faculty may want to use 
this project in their own classes, the problem used in class 
follows. 

Later this semester, we're planning to go to Dream 
Acres Farm, near Wykoff, MN, for a seminar on black- 
smithing. Although the day will likely be fun, going to 
school is not all about having fun - as the two semesters 
have progressed I hope that you've begun to realize that 
what we study in class and in lab is closely related to the 
phenomena we see in our lives every day. While there 
are many things we can think about in the context of our 
course material while at the farm, what I'd like for all of 
you to study is the efficiency of the heating process within 
the forge. 

To work and shape iron, the metal needs to be heated 
to temperatures well above the temperature at which green 
wood combusts. Iron is tempered, (based on the anneal- 
ing temperature the metal can have varying degrees of 
outer hardness, and inner flexibility) at a temperature 
range of 100 — 300°C, and forged (malleable because of 
the temperature and workable with tools) at temperatures 
of 700 — 1400°C. Because green wood starts to smoke 
at about 150° C , it isn't suitable for heating iron to these 
high temperatures. In fact, if you have a hot bed of coals 
and throw on a fresh chunk of wood, you'll actually cool 
the fire until the wood burns down to charcoal. Instead 
of fresh wood, we need to use a hotter burning fuel, like 
charcoal, coal, or coke (which is coal with the tar driven 
off by heat). There are a number of semi- scientific works 
available on blacksmithing , see for example, [26]. 

A piece of iron is heated by simply sticking it in a pile 
of fuel in the forge, and then forcing air through the fuel 
to increase the rate of fuel combustion. When the iron 
reaches the desired temperature, you remove it, work it, 
and then, when too cold to work anymore, reheat the piece 
in the forge again. 

The initial question is fairly simple. When a piece of 
metal is put into the fire to be heated, how efficient is the 
transfer of combustion energy into the metal? 

More specifically, please tackle the three following 
questions: 

1. If you take an automotive leaf spring, dimensions 
lj inches by -r inch by 12 inches, initially at 10°C, 
and heat it in the forge until the bottom 4 inches 
glow orange or hotter, about how much energy needs 
to be added to the iron for this to happen? 



2. If this heating takes about 45 seconds, what is 
the power input (per unit area) from the coals in 
the forge? What minimum amount of charcoal is 
burned to accomplish this? 

3. Use this estimate of energy flow rate in a numerical 
simulation of energy flow within the bar, and pre- 
dict a surface temperature distribution for the iron 
after 15, 30, and 45 seconds of heating. 

A sample solution to the first two questions is provided 
in Appendix B, and the theory underlying a numerical 
heat transfer scheme, which the students used to build a 
model, is provided in Appendix C. 

B. Numerical Preparation 

There are probably a number of ways to address this 
problem, it might be easiest at this point to stick a rod of 
iron (of known dimension and initial temperature) into 
the coals, heat the iron, and then pull it out and check 
the temperature profile with our digital thermometer. We 
should be able to figure out, based on the temperature gra- 
dient, how much energy went into the iron, and further, 
if we measure the change in mass of coal over the forging 
period, have a reasonable estimate for how much energy 
was given up in combustion. 

Before we actually go to the farm, I want you to build a 
computational model of the rod, heating in a pile of coals. 
When we actually go to the farm, while one person heats 
the rod, a second person can measure the time necessary 
for the metal to heat up to forging temperature. A good 
computational simulation of the process should be able to 
duplicate these results with reasonable accuracy. 

When we get back from the farm, I'd like you to refine 
your computational model to account for any discrepan- 
cies between the data you took, and the temperature pro- 
file you measured experimentally. 

C. While at the Farm 

While in the forge area, and after you've been in- 
structed about safety and technique, I want you to take 
a piece of metal stock and score regular marks on it with 
a file. Then heat up the piece of iron until its quite hot 
(red-orange is sufficient) noting the clock-time necessary 
for the piece to heat to this color. 

When heated, remove the iron from the coals and take 
temperature measurements with the IR thermometer at 
the marked locations, at 30 or 60 second intervals (what- 
ever seems feasible), for at least 10 minutes. 

Your measurements right after the iron is removed 
from the fire should correspond closely to the numerical 
heat transfer simulation you wrote before going on the 
lab. This experimental data set will allow you to calculate 
the rate at which heat flows from the coals into the fire 
(something we had to guess at in the simulation). 
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FIG. 1: (Color online) This image shows the basic set up 
of the lab. One benefit of the high temperatures involved is 
that a magnet can be used as a primitive thermometer. Iron 
undergoes a ferromagnetic/paramagnetic transition at about 
770° C, so by running a magnet down the length of the stock 
and finding this transition point, students can estimate how 
much heat has been transferred into the bar. For departments 
without an infrared thermometer, finding this magnetic tran- 
sition point is a straightforward (and cheap) way to obtain 
quantitative data to check a simulation. 

D. Reflection and Analysis after we return 

Given the data you took for the heating and cooling 
of the iron in the forge, refine the computational model 
for heat transfer you wrote before the trip. As a write up 
for the programming project please answer the following 
questions (including figures, equations, and diagrams as 
appropriate) : 

1. What average power density does a stoked fire pro- 
vide to a piece of iron stock? How did you figure 
this out? How accurate do you think this estimate 
is? (Include a well-reasoned estimate of the uncer- 
tainty in your answer.) 

2. Did you include radiative losses from the un-heated 
end of the iron in your model of heat transfer? 
Do you think they're necessary? Do you think the 
physics included in the computational model is com- 
plete? As always, justify any statements you make 
with respectable scientific arguments. 

IV. STUDENT SOLUTION 

The field trip ran roughly as described, and in addi- 
tion to the temperature measurements described, stu- 
dents were able to see the ferromagnetic/paramagnetic 
phase transition of iron by heating the bar and then 
bringing a large "cow" magnet near the sample, see fig- 
ure 1 . We also brought a small spectrometer along to the 




FIG. 2: (Color online) Students created a distance scale on 
the iron rods by cutting notches into bar stock with a file. 
There was a metal tape-measure handy in the blacksmith shop 
and students felt that inches were a natural distance unit to 
use in the lab. The simulations used metric units exclusively, 
and simulation output was converted into inches to make a 
comparison with the measured temperature data. 



shop and above the blackbody background, we saw clear 
sodium spectral emission lines at about 580nm. 

The remainder of this section (IV) is the solution pro- 
duced by one student, Nicole Schoolmeesters, who went 
on the trip. The code Nicole wrote for the project is 
available online [28]. She writes: 

At a blacksmithing workshop at Dream Acres, we in- 
vestigated heating iron rods to see first hand the differ- 
ent properties of iron. For this lab, we worked with an 
iron rod with the dimensions of 12 inches by l| inch 
by j inch that is initially at 10° C . A narrow notch was 
scratched across the width of the rod at one inch intervals 
to reference surface temperature measurements, see fig- 
ure 2. For each trial, approximately 4 inches of the rod 
was stuck into the hottest part of the coals for certain 
length of time. Then we quickly measured the tempera- 
ture with an IR thermometer at every inch notch mark 
along the length of the rod. Between each trial, the rod 
was quenched so that the initial temperature of it would 
return to approximately 10°C. The data recorded at the 
forge was then compared to the output data from the 
model rod in the simulation written before the field trip. 

The goal of this lab was to see if we could make our 
VPython program output temperature values along the 
surface of the iron rod that duplicates the real life data. 
The computer code simulated 30 seconds of heating, and 
the code's output could then be compared to the 15 sec- 
ond and a 30 second temperature measurements which we 
took at the blacksmithing forge. My program simulates 
the heat flow within a rod with the heat equation, which 
is derived from the ideas of calorimetry and heat conduc- 
tion, and is contained in Appendix C. When the rod is 
stuck into the inside of the forge, energy flows into the 
rod and is seen as an increase of temperature. I described 
this influx of heat as P CO ais > which is the the rate of en- 
ergy flow per unit area on the outside surface of the rod. 
For the portion not stuck into the coals, I also included 
in my simulation the non-contact heat flows from the rod 
to the rod's surroundings via the Stefan-Boltzmann ra- 
diation law (for a blackbody, emissivity 1). 

Although the simulation code is written with metric 
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FIG. 3: (Color online) Simulation and experimental surface 
temperature data for the steel (iron) rod shown in earlier fig- 
ures. The data was taken after the rod has been heated in 
the forge for 15 seconds. The purple, black, and red lines 
correspond to the simulation where heat is transferred into 
the rod from direct contact with the coals. The green line is 
from a modified simulation which adds an approximate radia- 
tion transfer scheme. This revised version of the code fits the 
experimental data (blue boxes) fairly well. Surface tempera- 
ture measurements on the rod were taken with a Fluke 68 IR 
Thermometer, which the manufacturer clams has an accuracy 
of±l°C. 



20 chunks along the l\ inch axis. This means that each 
"finite clement" in our model has dimensions of | inch 
by 32 inch by jg inch. 

Although the material properties of iron change with 
temperature, in this model we held the following prop- 
density, 



kg k- 



erties constant: specific heat, c p — 449 

P = 7870-^r; and thermal conductivity, k = 80.2-^. 
The entire system's temperature was initialized to 10°. 

The timestep for the program is determined by the nu- 
merical stability of the algorithm described in Appendix 
C, and was determined largely by trial and error. The 
largest timestep that we used is dt — 0.002 seconds. On 
a 2.16GHz MacBook the program then takes about 3 
hours to simulate 30 seconds of forge heating, [27]. 

The bulk of the code is derived in Appendix C, and an 
implementation in VPython is available as a supplemen- 
tary file on arxiv.org, [28]. 

A comparison of the forge temperature data and the 
prediction values from the simulation written before the 
trip is shown in figures 3 and 4. The calculated value 
for P coa is in the first trial was w 207 kW/m 2 . As the 
figures show, the simulation output and the actual forge 



FIG. 4: (Color online) This figure contains the same plotted 
variables as figure 3, but the time corresponds to 30 seconds 
of heating in the blacksmith forge. Note again the poor agree- 
ment of simulation to data when the radiation transfer scheme 
is not included. Figures 3-6 use a horizontal length scale, 
along the metal bar, of inches, where corresponds to the 
end of the rod in the coals. The simulation uses metric units, 
and output was converted to the english length in order to 
correspond to the measurements taken at the blacksmithing 
workshop. 



units, we used english units throughout the paper be- 
cause there was a carpenter's tape measure in the black- 
smith shop and inches was a natural unit for the prob- 
lem. In the simulation, the rod is subdivided into small 
chunks to more accurately model how energy flows into, 
through, and out of the rod. There are 48 chunks along 
the 12 inch axis, 16 chunks along the \ inch axis, and ^< 



data have a similar shape for the four inches of the rod 
stuck into the coals, but the simulation consistently pro- 
duces higher temperature values for both the 15 second 
intervals and the 30 second intervals. 

To understand why the simulation produced higher 
temperatures, I ran the simulation again and decreased 
the value of P, 



lis 



coals- The 15 second forge data suggested 
155kW/m 2 . The subsequent run has 15 second 



output which is better correlated to the 15 second forge 
data, but the simulation output was still too high for the 
30 second data (again, see figures Figures 3 and 4). Ac- 
cordingly, I decreased the heat input again, only this time 
using the 30 second forge data to determine a value of 
Pcoais ~ H5kW/m 2 . This time, the temperature output 
from VPython was too low for the 15 second line in that 
simulation, but the 30 second line shares a similar shape 
and is within the expected temperature range when com- 
pared to the forge data. From this set of varied inputs 
it seems that the rate of energy flow into the rod must 
decrease as the temperature increases within the rod. 

Although the 15 second forge data and the VPython 
output where P CO ais ~ l55kW/m 2 match well, another 
refinement to the VPython program could be made for 
the middle portion of the rod. When plotting the temper- 
ature data measured at the blacksmithing forge, there is a 
gradual decrease in temperature between the four inches 
of the rod stuck into the hottest part of the forge and the 
four inches at the opposite end of the rod that is unaf- 
fected by the heat from the forge. However, the VPython 
code was written so that there is an abrupt change to no 
additional heat inflow past the four inch mark where the 
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FIG. 5: The original model of thermal energy flow into an 
iron rod placed approximately 4 inches into a blacksmithing 
forge has a constant incoming power density, P coa is- The 
modified model decreases P coa is by 25% for every inch within 
the middle portion of the rod, as seen by the gradation of grey 
between the immersed 4 inches and the cold, outer portion of 
the rod. In both cases, the portion of the rod not inserted 
into the coals radiates heat to the colder surroundings via the 
Stefan-Boltzman Radiation law. 



rod is immersed into the forge. The middle portion of 
the rod, where it is still surrounded by hot coals but not 
physically touching the coals should have a decreasing 
amount of energy flow into it, proportional to the dis- 
tance a segment of the rod is from the surface of the 
bed of coals. This observation leads to the heat transfer 
scheme shown in figure 5. 

As implemented, this modification decreases P CO ais by 
25% for every inch within this middle area. For example, 
between inch 4 and inch 5, the power density is 75% of 
the original value. Between inch 5 and inch 6, the power 
density is half that of the original, and between inch 6 
and inch 7 it is 25% that of the original. The remaining 
part of the rod beyond inch 7 no longer has any heat flow 
into it from the coals. After including this modification, 
using P CO ais = 155kW/m 2 , as figure 3 shows, I found the 
the simulation output fit the forge data really well for the 
15 second temperature output compared to the previous 
trials. As for the 30 second data, shown in figure 4, the 
temperature output from the simulation was high for the 
first four inches, but fits the rest of the forge data fairly 
well. 

The thermal conductivity(fc), specific heat (c p ), and 
density (p) of iron change with temperature. Addition- 
ally, as a piece of iron is worked, the surface composition 
can change subtlety. Given this, I wondered if a different 
value of k — — , the coupling constant in the heat trans- 
fer equation, might also account for the discrepancy seen 
in the comparison of experimental data and the original 
simulation. (Basically, I wanted to see if the radiation 
transfer scheme described above was really necessary, or 
if variations in k might lead to the same surface tempera- 
ture distribution.) To test this idea, I ran the simulation 
again with a drastically increased and decreased value for 
k to see if there was a significant change to the shape of 
the graphs. As seen in figure 6, subtle variations in k 
cannot account for the differences seen between the tem- 
peratures measured at the forge and the temperatures 
produced by the simulation. 
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FIG. 6: To see if variations in specific heat capacity, den- 
sity, or thermal conductivity might account for the discrep- 
ancy between the original simulation and the experimental 
data, we re-ran the simulation (with the radiation scheme in- 
cluded) with wildly divergent k values of the original Ko ~ 

2.27 x KT 5 ^, K = 1.5/sn, and k = 0.5/t - As the output 
indicates, variations in « will not significantly affect the sim- 
ulation output. 



V. ALTERNATIVE IMPLEMENTATIONS 

While few campuses have blacksmith shops in their 
science complex [29], the direct comparison of simula- 
tion with experimental data is a powerful and necessary 
thing for an undergraduate in science to see and experi- 
ence. The data taken in this project is used to check the 
accuracy of the simulation and hopefully suggest other 
physics. Many other variations of this lab, not requiring 
a blacksmith's forge, are feasible and we now discuss a 
few options. 

The simulation described in Appendix C accounts for 
heat transfer within a box-shaped object with known 
power flux through the surface of the object. To simulate 
the forge, an instructor could use somthing a simple as 
an electric toaster oven or a "Weber" style charcoal grill, 
both of which would have power densities sufficient to 
heat a piece of metal stock. 

Although not discussed, the simulation should also 
model quenching of a hot metal object fairly well (al- 
though we have not taken data to verify this). Students 
could pre-heat a piece of iron in a water bath and then 
quench the stock in ice water or liquid nitrogen. In 
both cases, the fundamental physics should be similar, 
although a temperature dependent heat flux may be nec- 
essary if the object is quenched for a long time. 

As discussed in the caption to figure 1, iron undergoes 
a paramagnetic/ferromagnetic phase transition at about 
770°C, so in lieu of an IR thermometer, a simple bar 
magnet can be used as a cheap temperature probe by 
locating the point at which magnetism "goes away" and 
then tracking the position of this transition point as the 
bar is heated for different time intervals. Alternatively, 
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students can use digital cameras to record the color dis- 
tribution across a heated bar (as illustrated in figure 1), 
and then compare these images with the forging color 
scales availabe in blacksmithing books [30]. Of course, 
these color profiles are directly related to the blackbody 
distibution, and the opportunity for another learning cy- 
cle is again available. 

Finally, the magnetic phase transition described above 
makes a compelling follow-up project for interested stu- 
dents. Students can study the boltzmann factor (see [22], 
chapter 11 for a excellent intro-level derivation), and then 
create an Ising/Monte-Carlo style model of the iron stock 
which would give a qualitative picture of the phase tran- 
sition described above. 



VI. CONCLUSIONS 

Computational physics is a natural part of a learning 
cycle in physics, not only because it illustrates concepts 
which are hard to imagine, but also because comparisons 
of simulation to data can suggest (or demand) new phe- 
nomena. This paper describes such a problem, in which 
radiative transfer, initially neglected, needed to be in- 
cluded for the simulation to better model reality. The 
idea that data demands a better description of nature is 
the essence both of physics and inquiry. 

This problem could be adapted to other campuses by 
using a different heat source or sink in place of the black- 
smithing forge. For example, the power from the forge 
could be simulated by placing metal stock in a toaster 
oven. Similarly, the quenching process could be simu- 
lated by dipping stock, initially at room temperature, in 
a bath of ice- water or liquid nitrogen. 

In our view, the implementation details are less impor- 
tant than the feedback loop described explicitly in this 
work. The greatest beauty of physics is that it is funda- 
mentally a real description of the world. When creating 
computational models, students need to be reminded and 
persuaded of this connection to reality by evaluating their 
solutions with data sets which are as tangible as possible. 
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APPENDIX A: GRADING RUBRIC FOR 
PROGRAMMING PROBLEMS 

Rubric for Assessing student understanding of 
basic physics (Using computer modeling) 

This is the rubric that we will use to evaluate your 
computer programs. Of course, at a minimum your pro- 
gram must work, but to get a good grade for the assign- 
ment, you must adequately document your program with 
comments that explain what each part of the program is 
doing, explain why you are solving the problem in the 
way you chose, and explain the relevant physics . Some 
of these comments should be in the body of the code, 
and some of the documentation (like plots of the results, 
derivations, etc) should be in an attached report. As 
reflected in this rubric, you will be graded for not only 
how well your program works but how well you have ar- 
ticulated your solution/program/physical-model to the 
grader. 

APPENDIX B: HEAT FLOW ESTIMATES 

Question: If you take a leaf spring, dimensions li 
inches by 4 inch by 12 inches, initially at 10°C, and heat 
it in the forge until the bottom 4 inches glow orange or 
hotter, about how much energy needs to be added to the 
iron for this to happen? 

Solution: The calorimetry equation, dU = mcdT, can 
be used to determine the amount of energy needed to 
heat the first four inches of the rod. The mass is found 
by multiplying the volume of the rod that was heated 
(4 x 1 1 x j) by the density of iron. The approximate 
temperature of an iron rod that glows orange is 1000° C, 
so the change in temperature that the rod undergoes is 
about 1000°C. The amount of energy that is then needed 
to heat up the rod is the product of mass, specific heat, 
and the temperature change, about 74fcJ. 

Question: If this heating takes about 45 seconds, 
what is the power input (per unit area) from the coals in 
the forge? What minimum amount of charcoal is burned 
to accomplish this? 

Solution: If heated for 45 seconds, the power input 
from the coals can be calculated with the relation, ^H- = 
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Grading Category \ 
Score 


Below Standard/non- 
existent 
(0 = F) 


Standard 
(2 = C) 


Distinguished (4 = A) 


Translation of 
problem into a 
physical model 


Assumptions of the 
model are not 
identified 


Adequately identifies 
the assumptions 
made when making 
the model 


"Limits of 
Simulation's 

Correctly describe 
where you expect 
your model to not 
accurately reflect a 
real situation and 
what that might 
imply about your 
model. 


Description of physics 
principles used 


The physics principles 
identified. 


A correct description 
of the physics 
principles that you 
used. 


A correct description 
of why you used the 
principles that you 
used in the program. 


Algorithm used to 
calculate the results 
of your model. 


The numerical 
algorithm doesn't 
simulate reality 


The 

numerical algorithm 
is correct 


The 

numerical algorithm 
is correct, 
numerically stable, 
and computationally 
efficient 


Analysis of results 
(Evaluate the 
answer) 


The results from the 
program are not 
checked with physics 
principles. 


Uses the main 
principled) from the 
program to correctly 
and conceptually 
describe the results 
that the 
program should 
be calculating. In 
addition, you check 
your answer by 
looking at the 
extremes 


In addition to the 
main principle(s) 
describing the 
results, any new 
physical phenomena 
and/or interesting 
results are explained 
analytically in the 
report. 



FIG. 7: This rubric was developed by Andrew Ferstl and 
Nathan Moore when they started including VPython pro- 
gramming problems in their introductory physics classes. 
Rather than simply using a computer to generate a number, 
we emphasized in class the importance of trying to recognize 
and understand the physics that a computer model illustrates. 
In order to change the way that our students think about pro- 
gramming assignments, we gave them the following categories 
to guide their perspective, code, and solutions. 

Pcoais x Surface Area. This works out to be P CO ais ~ 
207^. 

The energy density of coke is about 29. 6M J /kg [31] 
and the energy needed to heat up the rod is 74fc J. There- 
fore, if all of the energy from the coke goes into the iron 
you'd need about 0.3kg of coke. Of course, this is a 
dramatic underestimate. 
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APPENDIX C: CREATING A NUMERICAL 
MODEL 

This appendix contains a derivation of the numerical 
method that students used to create a computational 
model for heat flow within a cubical object. 

To build a useful model, we need to understand the 
system fairly well. Lets assume that we're trying to build 
a wood chisel out of a leaf spring from a 1970 's era full- 
size pickup truck. This iron, (technically, steel), is high 
carbon, perhaps as high as 1.5%, and great for making 
tools. The piece of metal is roughly 2 inches wide, 16 
inches long, and about 1/4 inch thick. To make the stock 
malleable, we need to heat the bottom 4 inches of the 
spring by burying it in the pile of coke (or charcoal) in 
the forge. The question we want to answer is "If one end 
of the bar is stuck into the coals, after a 15 or 20 second 
heating period, what temperature distribution should we 
see along the length of the bar?" 



of the same material is called calorimetry. Formally, the 
definition relies on "heat capacity" , C, which is the rate 
of change of internal energy, U, with respect to temper- 
ature T. If we talk about the change in temperature and 
internal energy of a specific mass m of material, subject 
to a constant (atmospheric) pressure, the rate of change 
is called "specific heat", c p . The useful relation is then, 

AU = mCpAT, (C3) 

where the change in internal energy AU causes a change 
in temperature AT. Most metals have a specific heat 
of approximately 3i? per mole, a beautiful result from 
the Einstein Solid model. For iron, the specific value is 
c p — 449 k g K at room temperature, but there are signif- 
icant changes in this value as the the iron's temperature 
changes. 

2. How does the heat actually move in the solid? 



1. Necessary Basic Ideas 

We first need to remember a few basic relationships 
from thermodynamics to be able to explain the mech- 
anism for heat transfer within a material. The rate of 
heat flow, Q, through a barrier like the exterior wall of 
your house, is related to the temperature gradient be- 
tween the inside and outside of the house wall, AT, the 
area of the wall, A, and the thickness of the wall, Ax, by 
the equation, 



AQ 
At 



kAAT 
Ax 



(CI) 



In this equation, k is the "thermal conductivity" of 
the barrier. For air, k = 0.025 which is quite low, 
and should make sense, as still air does not conduct heat 
well. Iron has a conductivity of k w 50 — 80 J^ K [31], de- 
pending on temperature and composition. By contrast, 
copper has a k value of nearly 400 , which can also 
be explained in terms of common experience. 

Metal objects often feel cold because the high ther- 
mal conductivity of the material allows lots of heat to 
conduct away from your skin, which creates the sen- 
sation of "cold" . If you touch comparable sized pieces 
of iron and copper, it is likely that the copper will feel 
colder because if its greater thermal conductivity. On a 
hand-waving level, copper's greater conductivity occurs 
because there are more ways for energy to move within 
a metal, specifically electrical and acoustical conduction. 
Naturally, there is also a differential form of the heat con- 
duction law which is more suited to our purposes in this 
project, 



at ax 



(C2) 



The relationship between the increase in internal en- 
ergy of a material and the observed temperature change 



To model the heat flow, I suggest you use the classical 
heat equation, which can be derived by combining the 
previous results for calorimetry and conductivity. The 
derivation of this equation comes by using calculus to 
combine the formulas for energy storage (Calorimetry) 
and energy flow (Thermal Conduction). 

First, think about a small cubical chunk of material. 
The cube has density p, volume, V — AxAyAz, and 
mass, m = pAxAyAz. At a certain point in time (as- 
suming constant heat capacity), the change in internal 
energy of the box can be related to the change in the 
temperature of the box by the differential form of equa- 
tion C3 above, 



dU = c p pVdT 



(C4) 



In general, if the heat flow into or out of the box is 
described as Q, we can go a step further and say, 



dQ dU dT 

— : — = — = CpV — 

dt dt at 



(C5) 



Now, we can specify the heat flow into the box by 
thinking about the heat conduction equation described 
above. Specifically, if we think about ONLY the x direc- 
tion, and say that the box is centered on x, then there 



will be two relevant heat flows, one from the x - 



Ax 



side, 



and another from the x+ =|p side. Mathematically then, 
the net heat flow is the difference between what flows in 
and what flows out: 



dU 
dt 



= Q\x- 



(C6) 



We know the formula for heat conduction, equation C2, 
so the task now is plug in that relation and see the math 
that falls out. For the geometry we're using, A = AyAz 
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dU 
dt 



= kA 



dT 

dx~ L+£ 



dT 
~dx~^ 



(C7) 



Then, if wc combine the two equations for 4£r, C2 and 



C7, there's a beautiful result: 



dt 



dT 
~dt 



kA 
c p pV 



dT 

dx~ L+£ 



dT 
dx 1 



(C8) 



If you think about the definitions, 4 = a A a A a = 

J ' V AxAyAz Ax ' 

and then, using the fundamental definition of a deriva- 
tive, 



dT 
Itt 



k d 2 T 



c p p 



dx 2 



(C9) 



In three dimensions the arguments are the same (re- 
place x with y, etc), and the constant = n, and the 
result becomes: 



dT 

~dt 



f d 2 T 
\ dx 2 



d 2 T 
dy 2 



8 2 T 

dz 2 



(CIO) 



This final result is a "partial differential equation" , 
which can be solved analytically (ie using math to figure 
out a closed-form solution) for very simple geometries, 
and numerically for more complex (real) geometries and 
heat flows. It is a virtual certainty that many (or nearly 
all) of the appliances you plug into the wall are thermally 
simulated before being actually built. Obvious examples 
are CPU heat sinks, the cooling coils in refrigerators, the 
heat exchanger in a furnace, the engine bay of your car, 
etc. 



3. A Numerical Model of the Bar 

How do we use a computer to solve this problem? First, 
we'll have to think about some sort of array- type repre- 
sentation of the iron bar. One of the common approaches 
is to think about the material as consisting of little pieces, 
each with an individual temperature, conductivity, mass, 
etc. Given the size of Avagadro's number, we can't make 
the chunks of iron the size of an atom, or even a few 
atoms, but if the chunks are much smaller than the ob- 
ject itself, the results can be quite accurate. More specif- 
ically, we want the chunks of iron to be smaller than the 
spatial scale of the phenomena we're interested in study- 
ing. In this case, since the part of the bar plunged into 
the coals is about 4 inches long, we should use a chunk 
size that's much smaller than 4 inches. The obvious limit 
to this reasoning is that as the chunk size decreases, the 
computation time required for a simulation increases in 
proportion to the number of chunks. Lets assume that 



the bar has dimension L x x L y x L z , where L x = 16in, 
L y = l/4in, and L z = 2in. Further, lets say that all 
chunks have dimension Air x Ay x Az. You can imag- 
ine then that the number of chunks in the x direction is 
N x = L x /Ax and so on. 

With this discretization then, we can talk about a 
representation in computer code. Since the differential 
equation for heat flow is written in terms of the temper- 
ature, T, we can represent the continuous distribution 
of temperature as a 3-d array, temp[i] [j] [k] , where 
the indices run from i = 0, N x — 1, j = 0, N y — 1, and 
k = 0,Nz-l. 

Once the array of temperatures is defined, we can work 
with an individual array element by referencing it by its 
3-d location in the lattice of mass chunks. If we define the 
i = chunk to span from x — to x — dx, and the i — 
N x — 1 chunk to span from x — L x — dx to x — L x , then 
we can initialize the bar to have an initial temperature 
with an iteration over the whole array. The mechanics 
of this are specific to the computer language used, see 
the associated example program for an implementation 
in Visual Python. 

Remember, to use a computer to model the heat flow, 
we need to assume that the bar of iron is made up of little 
boxes, each with a temperature of temp[i] [j] [k] . The 
trick now is to re-write the equation for heat flow, equa- 
tion CIO, to accommodate this numerical approximation 
of the iron, rather than the continuous distribution that 
the equation above describes. 



4. Finite Differences 

This approach is inspired by a discussion of numeri- 
cal solutions to electrostatics in Chapter 5 of "Compu- 
tational Physics" by Giordano and Nakanishi[32]. Al- 
though the physical context is quite different, similar 
mathematical techniques can be used for both problems. 

To implement a computational solution, we need to fig- 
ure out how to represent in terms of temp [i] [j] [k] 
in the x, y, and z directions. To do this, you need to think 
about approximating this second derivative as a small 
box to box change, or "finite difference" . Specifically, if 
we want to know the spatial derivative of temperature, 
we can say that, 



dT 
dx 



AT 
Ax 



T(a 



Ax)-T(x) 
Ax 



(Gil) 



This approximation increases in accuracy as Ax gets 
smaller. With a little creative re-ordering of the x-axis, 
we can say that the derivative at x — x is, 



dT 

dx 



T(x + \Ax) - T(x ~ \Ax) 



(C12) 



x—xq Ax 

Now, if we want to describe the second derivative at 
x = xq, we just have to (creatively) say that the second 
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±\Ax. 



derivative is the difference between two first derivatives, 
and pick the first derivatives to occur at x = xq 
Mathematics should make this more clear: 

First, the second derivative, expressed in terms of first 
derivatives: 



d 2 T 



dx 2 



=x Ax V dx 



dT 



_ dT 

(C13) 

Next, the two first derivatives, evaluated at x — xq ± 



x— xq— h A a; 



\Ax 



dT 
dx 
dT 
dx 



_ T(x + Ax) - T(x ) 
i=i»+|Ai Ax 

^ T(x ) - T(x - Ax) 

x=x -hAx Ax 



(C14) 



If you combine the previous three equations, you'll have a 
nice approximation the second derivative of temperature 
in the x direction: 



d 2 T 
dx 2 



T(x + Ax) - 2T(x ) + T(x - Ax) 



(Ax) 2 

(C15) 

Now, if we want to move this approximation into python 
code, the trick is to remember that adjacent sites (i, i+1) 
are separated by dx or, Ax in the language of the deriva- 
tion. This means that the second derivative at the ith 
site can be represented as: 

ddT_dxdx = (temp [i+1] [j] [k] - 2 . 0*temp [i] [j] [k] 
+ temp[i-l] [j] [k])/(dx*dx) 

The other spatial derivatives can be similarly specified, 
see the example code for details. 



5. What about the edges? 

If you think about it carefully, you'll realize that the 
above formula for won't work at the x = or x = 
L x edges of the chunk of iron. Why? Well, there is no 
i+1 node if you're at the uppermost, i = N x — 1, value, 
and similarly, there is no i — 1 neighbor if you're at the 
i = location. To get around this problem, we can use a 
different approximation. Specifically, if we consider the 
simplest model of iron, in which the bar is completely 
isolated, there won't be a flow of heat to or from the 
outside, and so the flux would be zero across the edge. 
This means that the derivative normal to an edge is zero. 
The code that results (you can work it out from the above 
equations) is, 

# for the i direction: 
if (i==0) : 

ddT_dxdx = (temp [i+1] [j] [k] 



- temp[i] [j] [k])/(dx*dx) 
elif (i==(Nx-D) : 

ddT_dxdx = (-temp[i] [j] [k] 
+ temp[i-l] [j] [k])/(dx*dx) 
else: 

ddT_dxdx = (temp [i+1] [j] [k] 

- 2.0*temp[i] [j] [k] 

+ temp[i-l] [j] [k])/(dx*dx) 



6. What about the time derivative? 



As a reminder, the heat flow equation is: 



dT _ (d 2 T 



d 2 T 

dy 2 



d 2 T \ 
Ih 2 ) 



(C16) 



We've worked out a way to approximate the right side 
of this equation, which describes the spatial variation in 
the temperature (or internal energy) of the solid. The left 
side of the equation, describes the time rate of change 
of the temperature (or internal energy) . If we think about 
the movement of heat in the same way we've talked about 
the movement of particles in general mechanics, ic, if 
v = jjj, then x(t + At) = x(t) + vAt., then the simu- 
lation can progress as follows. Given an initial distribu- 
tion of temperature, temp [i] [j] [k] , once can compute a 
change to that distribution, d_temp[i] [j] [k] , and then 
say that the new energy distribution (after a time dt) is 
temp[i] [j] [k]+d_temp[i] [j] [k] . 



7. Putting it all together 

With everything we've talked about, from an initial 
distribution, temp[i] [j] [k] , you should be able to fig- 
ure out what the spatial changes are for every chunk of 
matter (ddT.dxdx, etc). Once these spatial changes are 
computed, you can figure out the change to each chunk's 
temperature (or internal energy) , and then finally, if this 
change is added to each chunk, a new distribution of en- 
ergy is obtained. With this machinery, you should be able 
to simulate the temperature throughout the material for 
as long a time as you'd like. 

A fancier (ie, more accurate) simulation would include 
effects like: the transfer of heat from the bar to the sur- 
roundings by radiation, the temperature-dependent heat 
capacity of iron, variations in the composition of iron (via 
the inclusion of carbon in the steel) along the bar, or the 
ability to simulate quenching in oil or water, etc. 



8. Heat flow from the forge, and from the rod 

The rate of energy flow per unit area, for something 
that's stuck into the forge is an unknown, and figuring 
out that value is one goal of this problem. If we call this 



energy flow (in units of -^), P coa is, then the effect on 
the system can be computed in the following way. 

If a certain mass clement (chunk) is on the k = sur- 
face of the bar, the surface area touching the outside, 
AxAy would permit the flow of energy Mathematically 
this means that the chunk of material has an additional 
d_temp[i] [j] [k] , which we can compute: 



— = P coals AxAy (C17) 

If we remember the definition of internal energy, AU — 
mcAT, we can solve for the additional temperature 
change, dT (or cLtemp), that occurs in the time inter- 
val At, because of the heat flow from the coals. 



dT 

pcAxAyAz — = P coals AxAy 
at 

dT Ax Ay At 

^£ * coals 



dt pepAxAyAz 

dT = +Pcoals ^L. (CIS) 
pepAz 

Finally, the iron rod will likely be warmer than the 
rest of the blacksmithing room, and so there is a second 
flux, this time from the rod to the surroundings, that 
is given by the Stefan-Boltzmann Radiation law (for a 
blackbody). According to this law, the power per unit 
area radiated by a hot object is given by, 



Pradiation — ®T (CI 9) 

where the power density is again in units of and 
where a = 5.6705 x 10~ 8 m 2^-4 - You should be able to 
use this relation to describe how heat leaves the surface 
of the rod which is not in contact with the coals. The 
specific implementation you could use (assuming again a 
mass element on the k=0 surface) is, 



dT — —Pradiation 7 — (C20) 

c p pAz 
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